Spontaneous symmetry breaking in active droplets provides a generic route to motility 
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We explore a generic mechanism whereby a droplet of active matter acquires motility by the 
spontaneous breakdown of a discrete symmetry. The model we study offers a simple representation 
of a "cell extract" comprising, e.g., a droplet of actomyosin solution. (Such extracts are used 
experimentally to model the cytoskeleton.) Actomyosin is an active gel whose polarity describes 
the mean sense of alignment of act in fibres. In the absence of polymerization and depolymerization 
processes ('treadmilling'), the gel's dynamics arises solely from the contractile motion of myosin 
motors; this should be unchanged when polarity is inverted. Our results suggest that motility can 
arise in the absence of treadmilling, by spontaneous symmetry breaking (SSB) of polarity inversion 
symmetry. Adapting our model to wall-bound cells in two dimensions, we find that as wall friction 
is reduced, treadmilling-induced motility falls but SSB-mediated motility rises. The latter might 
therefore be crucial in three dimensions where frictional forces are likely to be modest. At a supra- 
cellular level, the same generic mechanism can impart motility to aggregates of non-motile but 
active bacteria; we show that SSB in this (extensile) case leads generically to rotational as well as 
translational motion. 

PACS numbers: 



Living cells can move themselves around in a variety of 
different conditions and environments, and they exploit 
a range of strategies and mechanisms to do so. Uncov- 
ering the generic pathways to cell motility remains cen- 
tral to many important processes ranging from wound 
healing and tissue development [1 to immunological re- 
sponse and diseases such as cancer [2|. The best char- 
acterized case is that of a crawling cell on a planar 2D 
substrate or wall. Here motility is generally attributed 
to cytoskeletal actin filaments which polymerize at one 
end (+) and depolymerize at the other (— ) in a process 
called treadmilling. So long as the system has nonzero 
polarity P = (p) (where p is a unit tangent oriented from 
— to + and angle brackets denote an average over fila- 
ments), treadmilling leads to macroscopic motion. This 
exploits a Brownian ratchet mechanism in which forward 
fluctuations of the cell perimeter are locked in by polmer- 
ization [3] . However this mechanism requires a relatively 
solid anchor-point from which to propel the cell. For 
cells crawling on a 2D surface, this is provided by focal 
adhesions and other integrin rich structures [4j [5] . 

In vivo, cells often move in a three-dimensional envi- 
ronment such as an extracellular matrix and/or a tissue 
of cells [2j[6j[7]. Especially in quasi-spherical motile cells, 
both integrin-rich structures and mechanical anchoring 
are less in evidence, and the mechanisms of motility in 
such 3D environments remain unclear [8] . A recent study 
on tumour cells moving inside an elastic gel suggests that 
an important role in 3D locomotion may be played not 
by polymerization but by myosin contractility. This can 
lead to collective internal flows of actin that may ulti- 
mately propel the cell forward [7]. The contractile effect 
arises by a motor spanning two fibres causing an inward 
force pair (Fig. I]A) [9 j. This creates an active stress, 
usually modelled as cr^ lve = C,cP a P$ with c the concen- 
tration of active material, Greek suffices denote Carte- 
sian directions, and £ is an activity parameter (£ > 



for contractile systems). This raises an important issue 
of principle: how exactly does a tensorial active stress 
result in a vectorial propulsion velocity? 
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FIG. 1: (A) A minimal model for an active element or force 
dipole. The figure shows the quadrupolar nature of the hy- 
drodynamic flow around an extensile element or pusher (left) 
and a contractile element or puller (right). (B) Contractile 
stress is created when a motor protein (myosin) pulls protein 
filaments (actin) together in the cytoskeleton. (C) Polymer- 
ization of the actin filaments gives rise to an effective "self 
advection" velocity in the direction of the polarisation vector 
(orange arrow). 

Several studies use minimal models to address the 
fundamental physics of how activity imparts cell motil- 
ity. Experimental progress has focused on "cell ex- 
tracts" [lOj [n]: unregulated bags of cytoskeletal fila- 
ments (actin) and molecular motors (myosin), enclosed 
by a membrane. However, most of this work focuses 
on 2D crawling via the treadmill-ratchet mechanism de- 
scribed above. (3D systems are harder to study, and 
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selective inhibition of the treadmill dynamics is biochem- 
ically difficult [12 .) On the modelling side, generic the- 
ories have been proposed to make contact with the cell- 
extract data, again focusing mainly on 2D crawling and 
the treadmill-ratchet mechanism [T3HT5] . While the in- 
fluence of myosin contractility on cell shape during loco- 
motion has been addressed [14j [15] , in 2D this has not so 
far been thought sufficient by itself to lead to motility. 

Here we provide a detailed computational study of the 
effects of active stresses in a minimal 2D model of an acto- 
myosin cell extract. To confirm that our proposed motil- 
ity mechanism remains pertinent in 3D, we additionally 
perform selective (computationally intensive) simulations 
in that case. Our simulation model comprises a droplet 
of an active fluid or gel [T6l H7] , confined by interfacial 
tension <j, and surrounded by a Newtonian host fluid. 
Our model equations are based on established continuum 
precepts and outlined in Methods and Materials. Myosin 
contractility is represented by an active stress as detailed 
above. This term is invariant under global polarity inver- 
sion as are, with treadmilling absent, the full equations 
of motion (for further discussions of their symmetries, 
refer to Appendix). We show then that when the activ- 
ity parameter f exceeds a given threshold, an initially 
circular or spherical droplet spontaneously breaks that 
inversion symmetry, leading to an elastic splay of the po- 
larity field and to motion along ±P. This spontaneous 
symmetry breaking (SSB) manifests itself as a supercrit- 
ical Hopf bifurcation, which can alternatively be viewed 
as a continuous nonequilibrium phase transition. (The 
threshold value depend on both <j, and an effective elas- 
tic constant n penalising distortions of the ordered polar 
state.) In 2D, the resulting motile droplets have crescent- 
like shapes similar to some crawling cells [18] , whereas in 
3D our model predicts both spherical and concave shapes. 
Representing wall friction by a depth- averaged drag term 
in 2D, we find this slows down the S SB-mediated motion, 
but does not stop it altogether unless a critical drag is 
exceeded. 

We then introduce the treadmilling effect, which we 
represent by a self-advection parameter w. This breaks 
the global P <H> — P symmetry, but it does so directly, 
not spontaneously. (Note that in our model polarity is 
present even without treadmilling: see Appendix A for a 
discussion of the relevant physics.) This w is the speed 
at which, relative to the local suspending fluid (of veloc- 
ity v(r, £)), each filament is self-propelled along its own 
tangent. Even though asymmetric polymerization does 
not lead directly to mass transport, for actomyosin the 
resulting mass flux wcP should capture, in a highly sim- 
plified manner, the preferential growth of filaments by ad- 
dition of monomers at one end and loss at the other; see 
Fig.[lp. This simplified description is possible because we 
exclude the bath of monomers from the local mass den- 
sity c(r, t) of active material. (We assume that on aver- 
age, the monomers nonetheless keep up with the moving 
gel.) In general, self-advection leads to density gradients, 
which in turn cause hydrodynamic flows v(r, t) in the di- 



rection opposite to P(r,£). This backflow severely limits 
the effectiveness of self-advection in creating motility, al- 
most cancelling it for a droplet in free space as we discuss 
below. Our 2D study of the effect of wall friction shows 
however that high enough friction, by reducing v towards 
zero, restores treadmilling-induced motility with a speed 
that approaches w~P. 

Although our main focus is on actomyosin cell extracts, 
our framework provides a broader generic approach to 
droplet motility. We can thus investigate what happens 
when we reverse the sign of the active stresses to con- 
sider the extensile case, ( = — f > [HI [20]. The pri- 
mary experimental relevance of this case is to suspen- 
sions of bacteria, which push fluid out along their major 
axes and draw it in around the equator; see Fig. [T]A.. 
In this context, w is the bacterial swim speed; although 
this is nonzero for motile species, one may create vir- 
tually non- motile mutants (called "shakers") which still 
create extensile active stresses, e.g. by excessively in- 
creasing the tumbling rate [21]. A droplet of such or- 
ganisms can be created either by inducing an attraction 
between them (as our model effectively assumes) [22] or 
perhaps by confining them in an emulsion droplet. To 
attain nonzero P one further requires these particles to 
have net polar order (as opposed to a nematic state, for 
which there is orientational order, but equal numbers of 
particles with tangent ±p locally (see Appendix)). Set- 
ting aside the possible difficulties in meeting all those re- 
quirements experimentally, we predict that such "shaker" 
bacterial droplets could again break the symmetry and 
start moving spontaneously. Intriguingly, the predicted 
trajectories in this extensile case are more complicated 
than those of the contractile model. This is because, in 
extensile droplets, the S SB-mediated velocity forcing is 
in a direction perpendicular, rather than parallel, to the 
polarization vector P. Restoring nonzero w to describe 
the case of motile bacteria, the composition of these two 
motions leads to circular or spiralling trajectories of the 
droplet as a whole. 



I. RESULTS 

We first present results for the 2D contractile case, then 
briefly describe our findings in 3D, and finally give some 
further 2D results on extensile droplets. We initialized 
our simulation runs with a circular (or spherical) droplet 
within which the concentration of active material is taken 
to be a constant (c = Co), with c = outside. The polar- 
ization field P(r, t) within the droplet is initially uniform 
along the horizontal (x) axis and varies with concentra- 
tion as magnitude \ ( c ~ Ccr ) with c cr a threshold value 

V Ccr 

for polarization onset. We choose Co > c cr > so that 
|P| is nonzero inside the droplet and zero (isotropic) out- 
side, so that we assume the actomyosin network inside 
our active droplets is initially polarised. The system is 
then evolved via the equations of motion as specified in 
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Materials and Methods (discussed further in Appendix 
A), with chosen values of the activity parameter £ and 
the self- advect ion (treadmilling) parameter w. On a rela- 
tively short timescale both the internal concentration and 
the polarisation field relax towards the equilibrium values 
c eq and P eg (while still c = and P = externally) with 
some interfacial tension cr, set by minimization of our 
chosen free energy. (Our choice creates no anchoring of 
P at the surface so, in the absence of symmetry breaking, 
the polarization remains uniform; see Fig. [2]A. left. The 
effect of a soft anchoring is discussed in the Appendix 
A.) Having made one such relaxed droplet, the depen- 
dence of its behavior on £ and/or w was systematically 
explored by incrementing those quantities and waiting 
for steady state, before incrementing again. Because we 
are primarily interested in trends and symmetry break- 
ing phenomena, rather than quantitative predictions of 
where these will occur for specific materials, we report 
all results below in the natural units for lattice Boltz- 
mann simulations (LBU); the connection between these 
and physical units is discussed in Appendix D. 



A. Contractile stress can create motility via SSB 

Although contractile motor stress and actin tread- 
milling are generally both present in motile cells, it is 
illuminating to study these two mechanisms separately. 
We first consider a droplet with no external drag term 
(no wall friction) and no treadmilling term (w = 0), and 
vary the activity parameter C > 0. For low activity £, be- 
low some critical value f c , the droplet polarisation field 
P remains aligned uniformly along its initial direction 
x and the droplet remains stationary. However it be- 
comes slightly elongated in the direction perpendicular 
to the polarisation vector P as a result of the compe- 
tition between the contractile stress and the interfacial 
tension (see Fig. [2]A. middle). In this regime, the con- 
tractile stress set up a quadrupolar fluid flow around the 
droplet (Fig. [2^3 left), so that the whole droplet behaves 
as a large contractile element (compare Fig. [T]A. right). 
However it does not translate: there is no motility, and 
the droplet is a "shaker" rather than a "mover" [23] . 

As we increase f beyond £ c , the uniform polarization 
field P becomes unstable with respect to a splay defor- 
mation. This happens because the contractile stress is 
large enough to overcome the resistance to deformation 
mediated by the elastic constant n. The splay creates a 
state in which neighboring vectors P either fan outwards 
(V • P > 0) or inwards (V • P < 0). The first is shown in 
Fig. [2]A. right; the second is found by first taking its mir- 
ror image and then reversing P. This choice is made at 
random, spontaneously breaking the global polarity in- 
version symmetry. As soon as this happens, the droplet 
starts to move along the direction set by (V • P)P = ±P. 
This motion is attributable to the formation of a pair of 
flow vortices inside the droplet (Fig. [2^3 right). Such 
spontaneous propulsion is somewhat reminiscent of the 




0.001 0.002 0.003 0.004 0.005 0.006 
contractile activity £ 



FIG. 2: (A) Steady state configurations of a contractile ac- 
tive droplet without self- advect ion. The (red) arrows show 
the polarisation field P(r, t). Upon increasing the contrac- 
tile activity £, the droplet elongates perpendicular to P and 
then becomes unstable with respect to splay deformation at 
critical activity f c . When it splays, the droplet also sponta- 
neously moves in the direction of the green arrow. The time 
evolution of the system is shown in Supplementary Movie 1. 
(B) Left plot shows the velocity field of the droplet at f < ( c 
which is quadrupolar, like that around a contractile element 
(Fig. right). Right plot shows the velocity field of the 
splayed and moving active droplet which consists of two op- 
posing vortices. The boundary of the droplet itself is given 
by the dashed line. (C) Bifurcation diagram showing spon- 
taneous symmetry breaking from a uniform and stationary 
state to a splayed and moving state as the activity parameter 
f is increased. 



self-electrophoretic motion of a vesicle with active mem- 
brane pumps in an ionic solution [24j [25] . 

We note that in bulk active contractile fluids, the state 
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of uniform P is also generically unstable to splay fluctu- 
ations in ID [I9j HH1 H7] , which then lead to the onset of 
spontaneous flow. At one level, the S SB-induced motility 
transition described here can be viewed as a manifesta- 
tion of that bulk instability, albeit with two variations. 
First, spontaneous flows are present on both sides of our 
transition: as discussed above there is a quadrupolar flow 
field already in the non-moving state. Second, in bulk the 
critical activity level is nonzero only in finite systems, for 
which the transition is discontinuous, unlike ours (see be- 
low and Fig. (2p), and the resulting velocity field much 
more complicated [TTJ [20] . 

To illustrate our symmetry breaking motility mecha- 
nism more clearly, we plot the magnitude Vcm of the 
centre of mass velocity of the droplet V r cm as a function 
of activity ( in Fig. [2p. This bifurcation diagram shows 
a continuous non-equilibrium transition from a station- 
ary and uniform state to a moving and splayed state. 
Moreover, to within numerical accuracy the observations 
are consistent with a supercritical Hopf bifurcation, for 
which Vcm ~ (C — Cc) ' 5 - This mean- field like exponent 
is perhaps unsurprising as there is no noise in our simu- 
lations. Accordingly it might change in the presence of 
activity-generated noise [28], depending on whether the 
bifurcation remains low-dimensional or acquires a many- 
body critical character. 



B. SSB-induced motility is diminished by friction 

The results above are for a droplet in 2D surrounded 
by Newtonian fluid. To better describe experiments in- 
volving cell-crawling on a substrate, we now consider an 
additional frictional force between the solid wall and the 
cell. To do this within our 2D continuum model, we intro- 
duce an additional force density f 'friction = — 7 V to the 
momentum balance equation (see Materials and Meth- 
ods and Appendix). Here 7 is an effective coefficient 
of friction which depend on whether we have non-slip 
or partial-slip boundary conditions on the substrate and 
also on the thickness of the cell. It may represent con- 
ventional friction and/or a coarse-grained model of focal 
adhesions and other localized mechanical contacts. The 
presence of this friction will significantly quench the hy- 
drodynamic flow v. Since our SSB-induced motility re- 
quires a hydrodynamic vortex flow inside the droplet, the 
frictional force can greatly reduce the droplet's velocity, 
bringing it to rest beyond a critical level 7 C , at which 
point the symmetry is restored; the value of 7 C depends 
on activity and other model parameters. This is seen in 
a plot of Vcm against 7 in Fig. [3] 



C. Self-advective motility is enhanced by friction 

We now consider the case when there is only tread- 
milling (modeled as self-advection w > 0) but no con- 
tractile stress (£ = 0). As discussed previously, the 
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FIG. 3: Representative plots of droplet velocity Vcm against 
frictional coefficient 7 for motile droplets driven by: con- 
tractile stress only (solid red line) and polymerization/self- 
advection only (dashed blue line). The inset shows the 
polymerisation-only driven motility in the limit of large fric- 
tion in which the droplet velocity approaches "polymerisa- 
tion" speed wP. 



droplet will move along the polarization vector P with 
speed Vcm < w. Since there is no spontaneous sym- 
metry breaking involved, translational motion occurs for 
any nonzero w, in contrast to the threshold behavior seen 
for contractile SSB-motility. The most interesting aspect 
is the role of the friction parameter 7. Again we plot the 
velocity of the droplet Vcm as a function of 7 in Fig. [3] 
In contrast to the previous case, motility is significantly 
enhanced by the presence of friction. Indeed, in in the 
limit 7 — » 00, we have v = and Vcm — ^ (see Fig. [3] 
inset). 

The intersection of the two plots of droplet velocity 
versus friction (found respectively by switching off activ- 
ity or self-advection) defines a characteristic friction scale 
7. For 7 < 7 contractile stresses dominate cell motility, 
while for 7 > 7, self-advection is dominant. 

In most experiments on 2D crawling of cells/cell ex- 
tracts [5], the involvement of focal adhesions suggests 
that the high friction (treadmilling dominated) limit 
generically prevails. On the other hand, some recent ex- 
periments [29] directly identify spontaneous symmetry 
breaking of the actomyosin network as the initiator of 
polarized cell motility in keratocytes. Our work empha- 
sizes that spontaneous breaking of global polarity inver- 
sion symmetry arises from contractile motor activity, not 
from treadmilling. It is therefore arguable that the role 
of motor activity in 2D motility has so far been under- 
estimated. We note however that an equivalent discrete 
symmetry breaking would create motility if a pure tread- 
milling state of zero P a = (p a ), but finite nematic order 
((PaPp) — 7^ 0), spontaneously acquires polarity 

locally (see Appendix). 
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FIG. 4: Steady state conformations in 3D contractile droplets 
without self-advection on increasing activity f from (A) to 
(D). B and C are motile as indicated while A and D are 
stationary. (E) shows the toroidal fluid flow inside the motile 
droplet of steady state B. The time evolution of the droplets 
in (B) and (D) are shown in Supplementary Movies 2 and 3 
respectively. 



D. 3D droplets show a window of SSB-induced 
motility 



In the context of experiments on 3D tumour cells, it 
has been argued that motility is driven primarily by con- 
tractile stress [2j [7] , suggesting that the low-friction limit 
of our model prevails here. This accords with the much 
diminished part played by focal adhesions in 3D [8 j. For 
the 3D case we therefore neglect the friction term, and 
run selected simulations to confirm that the SSB route 
to motility remains operative. 

Fig. [4] shows steady state polarization fields inside a 3D 
contractile droplet with increasing values of the activity 
magnitude from A to D. As expected, the first steady 
state encountered is a symmetric but deformed immotile 
droplet (Fig.[4|\). The active stress contracts the droplet 
along ±P resulting in a lenticular shape. As we increase 
( beyond a critical £ c , splay instability spontaneously 
breaks symmetry and causes the droplet to move along 
(V.P)P just as in the 2D case. The droplet shape is con- 
cave and (as in the lenticular case) both it and the flow 
field resemble that created by rotating the 2D droplet 
about the P axis (Fig. [4^3). The resulting hydrodynamic 
flow (Fig. |4^1) therefore corresponds to a toroidal vortex 
ring. It would be interesting to see how these predictions 
compare with intracellular actin and fluid flow maps that 
might in future be measured for cells moving in 3D envi- 
ronments, for instance those studied in Ref. [2]. 

Interestingly, as ( is increased further, the droplet be- 
comes increasingly spherical (Fig. |4p), and finally sym- 
metry is restored, creating an immotile spherical droplet 
with a 'hedgehog' defect (of topological charge 1 as dic- 
tated by the polar ordering [30]) at the centre (Fig.[4p). 



E. Extensile SSB creates transverse or circular 
motility 



Our final results are for 2D extensile droplets (Fig.[T]A). 
Continuum descriptions of uniform extensile active flu- 
ids are widely used to describe dense bacterial suspen- 
sions pjJl [23] [31] . These results may therefore be rele- 
vant to bacterial droplets formed by aggregation in the 
presence of attractive forces [22], or possibly by confine- 
ment of bacteria within a droplet emulsion. Again, for 
fixed thermodynamic parameters controlling the elastic- 
ity and interfacial tension of the confined material, one 
can vary the activity parameter ( = — ( (now positive) 
and the self-advection parameter w which corresponds to 
the swimming speed of individual bacteria and is nonzero 
for "movers" but zero for "shakers" [19] [26] . 

Fig. [5}^ shows steady state configurations of a 2D 
purely extensile (w = 0) droplet at different ranges of 
activity. For ( < ( c i the droplet remains stationary but 
again elongates symmetrically, this time along the direc- 
tion of the polarization field ±P. The extensile stress 
creates the quadrupolar flow field expected of a large, 
extensile shaker. For ( beyond the critical value £ c i, the 
droplet again becomes unstable, but now with respect to 
bend deformation as opposed to splay. This gives rise to 
a horizontal vortex pair inside the droplet as opposed to 
the vertical one in the contractile case (Fig. [2^3 right). 
This flow field causes the droplet to move in a direction 
set by the sense of bending as P x ( V x P) (which can be 
upwards or downwards according to Fig. [5]A.). As can be 
seen from figure Fig. [5]A. this stationary-to-motile transi- 
tion can still be characterised as a continuous SSB transi- 
tion, however, the droplet speed Vcm attains a maximum 
at intermediate ( before falling to zero again. Beyond 
a second critical activity threshold (Cc2 in Fig. [5]A), the 
polarisation pattern oscillates continuously while remain- 
ing symmetric at all times and consequently, the droplet 
again becomes non-motile. 

Turning finally to the case of nonzero self-advection 
(w > 0) we find that this, combined with the extensile 
motility, can give rise to an intriguing spiralling motion 
(see Fig. [5^3). This arises because the SSB-induced motil- 
ity is at right angles to the polarization (Fig. [5^3 right). 
On the other hand, self-advection everywhere transports 
material along P(r, t) locally; since the polarization field 
is spontaneously curved, this by itself would give a cir- 
cular droplet orbit. Compounding these two motions 
typically leads to outward spiral trajectories as shown 
in Fig. [5j3. This outcome contrasts with the contractile 
case where both the SSB-induced and the self-advective 
motion (the latter averaged over the droplet configura- 
tion) point either together or oppositely along the ±P 
direction and only straight line motion can result. 
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FIG. 5: (A) Plot of centre of mass velocity against activity 
(£ = — C > 0) for extensile droplet without self-advection. It 
shows continuous transitions from stationary to motile and 
then from motile to oscillatory at critical activity £ c i and £ C 2 
respectively. Also shown are the steady state polarisation field 
P for the stationary and motile case. The movie of the time 
evolution of the system is shown in Supplementary Movie 4. 
(B) The presence of both extensile stress and self-advection 
leads to an outward spiral trajectory (solid green lines). Also 
shown are the snapshots of the polarisation field at different 
timesteps (red arrows). 



II. DISCUSSION AND CONCLUSIONS 

Our simulations of contractile droplets can be viewed 
as a simple in silico analogue of the in vitro cell-extract 
studies that have been used to dissect the biophysical 



ingredients of motility. Our simulation work powerfully 
complements these studies, by allowing us to isolate the 
role of contractile (motor) activity in cell locomotion: 
this is very difficult in the laboratory, where current 
strategies for inhibiting polymerization dynamics (tread- 
milling) can severely impair other key subcellular pro- 
cesses [T2] . 

Our approach likewise complements those of previous 
theories [14] [15] [32] which have concentrated on tread- 
milling as the main driver of motility. We have shown in 
2D that contractile stresses alone can not only shape the 
rear of a crawling cell (Fig. [2|, but also create motility 
itself, provided that the motor forces are large enough to 
create an asymmetric circulatory flow as in Fig. [2j3. In- 
terestingly, keratocyte cells crawling on glass seemingly 
do exploit myosin activity to set up an intracellular actin 
flow in the rear of a cell which ultimately polarises it and 
makes it motile [29]. Nonetheless, our study of the the 
effects of a frictional term, which promotes the motility 
created by treadmilling but inhibits that caused by con- 
tractile stress, lends support to the view that cell crawling 
on a wall is usually dominated by the treadmilling. 

In contrast, in 3D cell motility, recent work suggests 
that treadmilling plays at most a minor role |2j [7] . In 
the 3D case, therefore, our work describes a simple and 
compelling mechanism for how spontaneous translational 
motion can in principle arise solely by the action of a 
contractile stress. This requires spontaneous symmetry 
breaking, mediated in our case by splay deformation in 
response to that stress. Our 2D and 3D simulations go 
beyond the 2D theory of Ref. [7 by addressing the dy- 
namics of the polarization field. We do however make 
some important simplifications: our droplets are con- 
fined only by interfacial tension not by an elastic mem- 
brane; we treat treadmilling as a simple self-advection; 
and we do not address any direct transition between ne- 
matic and polar order, despite assuming polarity inver- 
sion symmetry at thermodynamic level. Improving the 
model in these respects will require a more detailed mi- 
croscopic derivation which we shall leave to future work. 
To test whether our model indeed captures the biophysics 
of 3D cell motility, it would be exciting to visualize ex- 
perimentally the detailed cytoskeletal organization and 
flow fields, e.g. for cells moving through matrigel [2j[7]. 

Our generic framework is not limited to contractile ac- 
tomyosin networks. Indeed we have discussed the case 
of extensile droplets, possibly relevant to aggregates or 
emulsions of active but immotile bacteria; here transla- 
tional motility arises by spontaneous symmetry breaking 
only at intermediate activity, and is mediated by bend 
rather than splay deformation. The addition of self- 
advection along the bent polarization field then leads in 
addition to rotational motion. We note that rotational 
and translational motility of small bacterial aggregates 
was recently observed, but attributed to a somewhat dif- 
ferent mechanism where symmetry is broken by frozen-in 
statistical fluctuations rather than SSB [22] . 

Finally, our hydrodynamic equations of motion, or 
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close variants of these, might in some cases be applicable 
to concentrated eukaryotic cell masses such as tissue [33] . 
Within a tissue each cell exerts forces on its neighbors 
which at the lowest order continuum level creates a cer- 
tain density of force dipoles [33] ; the velocity field v then 
describes the slow migration of cells inside the tissue. 
It is intriguing to note that the large-scale tissue flow 
in animal cells during gastrulation may break the sym- 
metry to form vortices, similarly to our active droplets. 
This is the case of the so-called "polonaise movements" 
which are observed in the developing chick embryo pQ, 
and which are important to form the correct supercel- 
lular structure. In this context a "polarisation" field is 
sometimes used to describe the orientation of individual 
cells |53]. The relation, if any, between the onset of this 
vortex flow and that seen in our droplets remains to be 
explored. 

III. MATERIALS AND METHODS 

We briefly outline here the hydrodynamic model used 
in this work (more details are in Appendix). We consider 
a fluid, comprising a mixture of active material and sol- 
vent, with constant total mass density p. The hydrody- 
namic variables whose dynamics we monitor are: (i) the 
concentration of active material c(r, t), (ii) the fluid ve- 
locity v(r, t) (with incompressibility requiring V.v = 0), 
and (hi) the polarization field P(r,£) = (p) as defined 
previously. 

Although an active droplet is a nonequilibrium system, 
we introduce the following free energy functional to de- 
scribe its equilibrium physics in the passive limit of zero 
activity: 

F[c,P] = J 'd*r{V{c)+ k -\Vc\> -f^l\P\\l) 

+ ||P| 4 + |( V P) 2 } 

Here a > is a phenomenological free energy amplitude, 
k determines the droplet interfacial tension, and n is an 
effective elastic constant. This choice of F[c, P] leads 
to a continuous isotropic-to-polar transition at c = c cr . 
To confine the active material into a droplet, we choose: 



^( c ) = 4^c 2 (c — Co) 2 and set cq > c cr . This creates 
two free energy minima corresponding to a phase of pure 
passive solvent (external to the droplet, c = and P = 0) 
and a polar active phase (inside the droplet, c = c eq > c cr 
andP = P ec? ). 

Treating the active material as locally conserved, the 
time evolution of the concentration field c(r, t) can then 
be written as a convective-diffusion equation: 

Be f SF\ 

^+V - (c(v + ^P)-MV— J =0 (2) 

where M is a thermodynamic mobility parameter and w 
is the self-advection parameter. 

The dynamics of the polarisation field P(r, t) follows 
an "active nematic" evolution [20], given by 

BY* 1 W 

_+(( V+w p).v)p = -e-p+^- p -fsp (3) 

where v and Q are the symmetric and anti-symmetric 
parts of the velocity gradient tensor Vv. T is the rota- 
tional viscosity and £ is related to the geometry of the 
active particles [15] . 

Force balance is ensured through the Navier- Stokes 
equation, p (d/dt + v • V) v = -VP + V • g total - 7V, 
where P is the isotropic pressure, —7V is the friction 
force per unit volume and a total is the total stress in the 
fluid which includes viscous, elastic/Ericksen, interfacial 
and "active" stresses (see Appendix for details and for a 
discussion of the effects of additional terms, allowed by 
symmetry, in the equations of motion). The active stress 
is (j a a f ve = (cP a Pp [23] where C > for contractile ac- 
tivity and ( = — ( > for extensile. 

To solve these equations in 2D and 3D, we performed 
hybrid lattice Boltzmann simulations, as done previously 
for other active flows [TTj [20] . 
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Appendix A: Hydrodynamic description 
1. Choice of order parameters 

We model actin filaments in the cytoskeleton, or bac- 
teria in suspension, as rod-shaped active polar particles. 
In the hydrodynamic (continuum) limit, the dynamics of 
these dense suspensions of active polar particles can be 
described by a few continuum variables. 

Our chosen hydrodynamic variables are: the concen- 
tration of the active particles c(r,t); the average velocity 
field v(r, t) of both the active particles and the solvent; 
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and finally the polarisation field P(r,t) which is defined 
as a mesoscopic average orientation of the polar parti- 
cles: P = (p) (with the average taken over molecular 
orientations p; see main text). The total mass density 
is assumed to be constant throughout (so the fluid is in- 
compressible) . 

Importantly, in this paper we assume P to be nonzero, 
even in the absence of self-advection (w = 0). We believe 
this to be the correct description for actomyosin gels in 
a cytoskeletal context: that is, in any mesoscopic region 
containing N aligned actin fibres, we assume an excess 
of one orientation over the other that scales as A 7 ", not as 
TV 1 / 2 (as would arise for random orientations) [35] . More- 
over, relaxation of P in the absence of activity requires 
slow reversal of filament directions; therefore the local 
polarization will not suddenly disappear if all activity is 
abruptly switched off. The resulting presence of a quasi- 
static polarization is important in principle, because it 
restricts the nature of the allowed orientational defects 
in the system to those of integer topological charge [30] . 

An alternative theory would suppose that without self- 
advection the system is generically in a nematic state, 
which has P = but a nonzero tensor order parame- 
ter Q = (pp — 1 1). This describes a different situation 

in which half-integer defects become possible (at least 
when w = 0). A brief discussion of the onset of sponta- 
neous motility in a nematic model is given in Appendix 
B. Meanwhile we emphasize that orientation of actin, as 
distinct from self-advection of any kind, is what the po- 
larization field P represents in our model. 

Note that in some parts of the liquid crystal litera- 
ture [36-38 , the term "nematic" is used for a system 
whose mathematical description involves a nonzero vec- 
torial order parameter (say P), but whose governing 
equations respect the global symmetry 



P(r)->-P(r) Vr 



(Al) 



(as will arise for f3 = w = in the equations to be devel- 
oped below). We believe this terminology in our specific 
problem is best avoided: it gives the impression that the 
spontaneous breakdown of this global symmetry is equiv- 
alent to a genuine transition from nematic (Q ^ 0,P = 0) 

to polar (P 0). But mathematically there is no such 
equivalence. Indeed, the nematic-to-polar transition is 
from a state in which half-integer defects are allowed, 
to one in which they are forbidden. In contrast, spon- 



taneous breakdown of the symmetry Eq. Al represents 
a transition between two states that both forbid half- 
integer defects. Put differently, the global symmetry in 
Eq. |A1| is not equivalent to a local symmetry, in which 
one can freely reverse the orientation p of all rods within 
any arbitrary subvolume. That local symmetry is the 
defining feature of the nematic state. It is fully built into 
any description based on the tensor order parameter Q, 

but the local symmetry is hidden if one replaces this with 
a director field n (a unit vector whose orientation is de- 
fined by the major principal axis of Q). For instance, a 



fictitious discontinuity (domain wall) between two slabs 
with n x — ±1 (say) correctly has zero energy cost in the 
tensorial description but acquires a spurious, and indeed 
infinite, square gradient energy ~ k J (\7n) 2 dV in the di- 
rector field approach. The corresponding domain wall 
would, in contrast, be physically real in a system with 
polar order. 



2. Hydrodynamic equations 

The time-evolution of hydrodynamic variables in ac- 
tive systems can either be derived analytically [19] or 
phenomenologically [23] - we will follow the latter av- 
enue in this work. 

In deriving phenomenologically the hydrodynamic 
equations, we first introduce the free energy functional: 



F[c,P] = J d 3 r{V(< 



+ ||Vc| 2 -f(^ 
-(VP) 2 + /3P- Vc}. 



|P#2) 



In the passive limit the system tends to minimise this 
free energy. Note that we have assumed the single elas- 
tic constant (= n) approximation. Furthermore, within 
our theory in the passive limit there is a second order 
transition from isotropic (|P| = 0) to polar (|P| > 0) 
phase at critical concentration c cr . Other phenomeno- 
logical parameters are a, which controls the isotropic-to- 
polar transition, and /c, which in conjunction with V(c) 
determines the interfacial tension. 

To enable droplet formation, we choose: 



V(c) 



4c! 



c 2 (c-c ) s 



(A3) 



and set cq > c cr . This creates two free energy minima 
corresponding to a phase of pure passive solvent (external 
to the droplet, c = and P = 0) and a polar active phase 



(inside the droplet, c 



> c cr and P = Peg)- Then 



c eq and P eq are found by minimizing F[c, P] in a state 
of uniform c and P. The resulting surface tension a will 
depend on a, /c, and the elastic constant k. In the special 
case of c cr = cq/2, the expression for c eq and P eg are: 



Co Co a 
+ TV 1+ 2^ 



eq 



1/4 



(A4) 
(A5) 



where p is a unit vector. 

Finally, the last term in Eq. |A2[ /3P • Vc, represents 
a soft anchoring of P to the droplet interface (so that P 
will tend to point outwards at the droplet perimeter for 
j3 > 0). For simplicity this was set to zero in the main 
text. Note that for given V(c) and fc, j3 determines a 
dimensionless ratio H = (71 — 72)/ (71 + 72) where 71,2 
are the interfacial tensions with polarity directed along 
the outward or inward normal. 
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A nonzero value of H (and hence of (3) is allowed by 
symmetry and indeed in a passive system of simple am- 
phiphilic molecules one would expect H value of order 
unity. (The interfacial tension is quite different for am- 
phiphiles correctly oriented at an interface than for ones 
pointing the wrong way.) However, for an actomyosin 
droplet, much smaller values of H appear likely: indeed 
we see no obvious mechanism to suggest a strong pref- 
erence of either the positive or negative end of an actin 
fibre for the droplet surface. Accordingly in the main 
text we set f3 = 0, so that the free energy is in varia nt un- 
der the global symmetry transformation of Eq .jAl| (This 
can still be broken by w.) However, in Fig. p] we show 
the effect of adding a small f3 term on the spontaneous 
symmetry breaking transition that leads to motility in 
the absence of self-advection (w = 0). Here, we choose 
f3 = 0.0001 which roughly corresponds to H ~ 10 -4 . The 
effect resembles that of applying a weak external field to 
a ferromagnet undergoing an Ising-type phase transition. 
Just as in that case, we can expect a detailed understand- 
ing of the zero-field (j3 = 0) case to offer fruitful mecha- 
nistic insights, even if in reality a small field is generically 
present. 

Since the total number of active particles is conserved, 
the time evolution of the concentration field c(r, t) can 
now be written as a convective-diffusion equation: 



dc 

at 



V- c(v + ^P)-MV 



SF 







(A6) 



where M is the mobility of the active particles, related 
to the diffusion constant by D ~ Ma and w\ is the self- 
advection or the speed of the active particles relative to 
the bulk fluid. In the main text, we stated that w\ (there 
denoted w) is related to the velocity of actin treadmilling 
(which is, in turn, proportional to the rate of actin poly- 
merisation at the positive end of the filament). However, 
a contribution to w\ may also arise due to other active 
processes in actomyosin systems such as when the motor 
proteins preferentially pull the actin filaments towards 
ther positive ends more favourably. Note also that the 
5F/Sc term is the chemical potential derived from the 



free energy in Eq. A2 



The dynamics of the polarisation field P(r, t) is bor- 
rowed from polar liquid crystal theory, and can be written 
as (see [41 for a more thorough derivation and justifi- 
cation) 



<9P 

~dt 



-((v + w 2 P)- V)P 



1 SF 

-g.p + &.p- f _ (A7) 



where v and Q are the symmetric and anti-symmetric 
parts of the velocity gradient tensor Vu. T is the ro- 
tational viscosity and £ is a shape factor related to the 
geometry of the active particles: £ > for rod-like par- 
ticles and £ < for oblate particles. Here we take £ 
positive as seems appropriate for filamentary contractile 
matter such as actin networks. The parameter £ also 
determines whether the particles are shear-aligning (for 



|£| > 1) or shear-tumbling (for |£| < 1). In this paper we 
assume £ > 1. Note that Eq. |A7| contains in principle a 
second self-advection parameter that, for an arbitrary 
model of activity, need not equal w\. However, it seems 
highly plausible that (at least for a treadmilling mecha- 
nism) both concentration and polarity should advect at 
the same rate; in this work we therefore set w\ = W2 = w. 
Note in addition that an even more complete approach 
would allow two furt her terms, i^3V(P 2 ) and k; 4 P(V.P) 
to appear in Eq. A7 [42 . Inclusion of these terms could 
lead to a significant rounding of the transition along the 
lines discussed above for the j3 term. 

Force balance in our system is enforced through the 
incompressible Navier-Stokes equation, 



V- v 







(A8) 



p [ — + v • V ) v = -VP + V • g totaL - 7 v (A9) 

where P is the isotropic pressure and —7V is the friction 
force, already discussed in the main text. Here a total is 
the total hydrodynamic stress which includes the active 
stress. There are four contributions to the hydrodynamic 
stress: 

^.total ^.viscous _|_ ^.elastic _|_ ^.interface _|_ ^.active ^^^Qj 

The first one is the viscous/dissipative stress which can 
be written as: a™p cous = r](d a v^ + dpv a ) where 77 is 
the shear viscosity and the greek indices indicate carte- 
sian coordinates. Next there is the elastic/Ericksen stress 
borrowed from the liquid crystal dynamics: 

<T e J$* tiC = ^(Pah-P^h a )-^(P a h^P^h a )-Kd a P^P y 

(All) 

The interfacial stress between the active phase and the 
passive phase is 



interface 
J a(3 



df 



d (d Q c 



-d a c (A12) 



similar to that of binary fluid [39 , with / defined to be 
the free energy density. The active stress can be derived 
by summing the contributions from each force dipole and 
coarse-graining [23] which results in: 



-active 
T a(3 



= (cP a Pp 



(A13) 



where f is the activity parameter which is positive for 
contractile particles (pullers) and negative for extensile 
particles (pushers). Note that active stress contributions 
proportional to d a Pp are also allowed by symmetry in 
principle [43 . These higher order terms are plausibly 
relevant to bacterial swimmers which may have a high 
degree of asymmetry between the particle's head and tail. 
However, for actomyosin solutions (comprising elongated 
fibres of uniform width, with the activity provided by di- 
lute motor proteins crawling along those fibres) any such 
terms are likely to be small and we omit them from our 
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model. The magnitude of the activity £ is proportional 
to the strength of the force dipoles. The active stress is 
fundamental in our theory, as it is this term which drives 
the system out of equilibrium. 




0.001 0.002 0.003 0.004 0.005 0.006 



contractile activity £ 



ours but which is much more closely related to it, is 
when a truly nematic active droplet (Q ^ 0, P = 0), 
spontaneously breaks spatial symmetry to create a splay 
or bend deformation. For a uniaxial nematic, Q a p = 
S (n a np — with n the director field; so long as de- 

fects are not involved, the dynamics in this case should 
be extremely similar to a polar system with P = S 1 ^ 2 ^ 
[30] . For instance, Fig. [7] shows one example of a spon- 
taneously moving state in a droplet of an active nematic. 
Here, the splay instability due to contractile stress creates 
a vortex flow inside the droplet which is indeed similar 
to the polar droplet case discussed in the main text. The 
droplet shape is somewhat different, but not dissimilar to 
the ones found for the polar droplet at somewhat higher 
activity values than shown in Fig. 2(B) of the main text, 
whose 3D analog is seen in Fig. 4(C). 

On the other hand, as shown already in Fig. 4(D) 
of the main text, defects can and do arise in our ac- 
tive droplet system. If, as we claimed above, nonzero P 
is the generic situation for actomyosin gels even under 
quasi-static conditions and with w = [40], then our 
description based on the dynamics of P is preferable on 
principle to one based on Q. 



(A) 



FIG. 6: Plots of the droplet velocity as a function of contrac- 
tile activity for /3 = and f3 = 0.0001 (roughly corresponds to 
H rsj 10 -4 ). This has a similar effect to a weak external field 
in para-/ferromagnetic transition which destroys the second 
order transition at the critical point. (Compare Figure 2(c) 
of main text.) 



Appendix B: Alternative models 
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As noted previously, in the absence of self-advection 
(w — 0) and with no anchoring term (/3 = 0), the equa- 
tions above are invariant under the global symmetry of 
Eq. |A1| This is spontaneously broken at large enough 
activity \C\ to create a new route to motility, as explored 
in the main text. 

As already explained, when defects are allowed for, 
the breaking of a global P — » — P symmetry (which 
leads to motility in our case) is mathematically distinct 
from a nematic-to-polar transition. Therefore one pos- 
sible but inequivalent route to motility, already men- 
tioned in the main text, is where a truly nematic system 
(Q / 0,P = 0) becomes polar (P / 0). If this happens 

in a system where treadmilling (say) of individual rods is 
already present at the molecular level, then this transi- 
tion should cause the onset of macroscopic self-advection 
with speed w oc |P|. This describes a mechanistically dif- 
ferent way in which spontaneous symmetry breaking can 
lead to motility; however no motility would ever arise, by 
this route, in the absence of self-advection (w = 0). 

A way in which spontaneous motility can arise even 
with w = 0, which remains mathematically distinct from 



FIG. 7: A droplet of active nematic (Q a p > and P a = 0) 
can also become motile due to the active stress: a a p = C c Qai3- 
(A) shows a typical configuration of the director field in a 
contractile droplet (£ > 0) moving to the right. (B) shows 
the corresponding velocity field consisting of a pair of vortices 
inside the droplet. The dashed line represents the droplet 
interface. 



Appendix C: Simple scaling analysis 

In the absence of self-advection, we observe a 
stationary-to-motile transition at some critical activity £ c 
which is accompanied by spontaneous symmetry break- 
As described in 



ing of the global symmetry in Eq. Al 



the main text, below the critical activity £ c , the droplet 
remains stationary while the polarisation field is uniform 
inside the droplet. In addition, the droplet is also elon- 
gated due the competition between the active stress and 
the interfacial tension. In general the shortest radius of 
the elongated droplet is a function of the activity £ and 
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surface tension a, or R((,a). 

At large enough activity, above the critical value £ c , 
the droplet becomes motile. This transition may be un- 
derstood as driven by spontaneous elastic deformations 
(bending or splay) - therefore the critical value ( c can be 
estimated by equating the active and elastic stresses, as 
follows 



Appendix D: Lattice Boltzmann units 



^.active ^.elastic 

a ~ cr 



(CI) 



Eqs. (All ) and ( A13) lead to g active ~ ( Co and g elastic 
hence to the following scaling law: 



Cc 



c R(( c ,a) 2 



(C2) 



which represents a mean-field estimate of the location of 
the critical point. Notice that the surface tension a will 
also depend on n in general. This scaling law also has the 
same form as that of ID spontaneous flow transition [44] . 
To test this scaling law, we may plot the droplet speed 
Vcm as a function of a dimensionless quantity: 



Cc R(C,a) 2 



(C3) 



near the critical point for different values of elastic con- 
stant n and initial radius of the droplet Rq = R(( = 0, a). 
These plots are given in Fig. U We can see from the fig- 
ure that all the four curves corresponding to different pa- 
rameters have almost the same critical point at (j) ~ 14.5, 
thereby validating our approximate scaling analysis. 



R ■ 15 and k = 0.04 
R = 16 and k = 0.04 
R = 17 and k = 0.04 
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FIG. 8: Plots of the droplet speed as a function of dimension- 
less parameter cj) for different values of initial droplet radius 
Ro = R(£ = 0,a) and elastic constant k. 



To establish an approximate correspondence between 
the natural simulation units (lattice Boltzmann units) 
and those of typical cell extract experiments, we choose 
the length-scale, time-scale and force-scale to be: L = 
1 /im, T = 10 ms, and F = 100 nN respectively. In lat- 
tice Boltzmann/simulation units these three scales are 
all defined to be equal to one exactly: that is, the lattice 
parameter is taken as the unit of length and the timestep 
the unit of time. The variables and parameters used in 
our model are given in the Table 1 with their respective 
values in simulation and physical units. These parame- 
ters have been chosen to be roughly in the same order 
of magnitude to the ones quoted in [5j |46]. Note that, 
in line with standard practice in Lattice Boltzmann sim- 
ulations, we chose the fluid mass density p to be much 
larger than the actual mass density of a real solvent (wa- 
ter) [45 . This is acceptable so long as the role of in- 
ertial (characterized by the Reynolds number) remains 
small; the procedure speeds up the computations by sev- 
eral orders of magnitude. In this case the choice of the 
force/density scale above gives the Reynolds number of 
about Re ~ 0.03 which is small enough for the flow to be 
laminar though much larger than the actual experimental 
value which is of the order of Re ~ 10 -12 . 
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TABLE I: * 

Table. 1. Typical values of the physical quantities used in the simulations. This choice of parameters are made to be 
consistent to other physical estimates in [5j H6] . 



Model variables and parameters 


Simulation units 


Physical units 


Effective shear viscosity, r\ 


5/3 

/ 


1.67kPa/s 


Effective elastic constant, n 


0.04 


4nN 


Shape factor, £ 


1.1 


dimensionless 


Self-ad vect ion/polymerisation speed, w 


0.0015 


0.15 jLtm/s 


Average concentration of actin- myosin pairs, Co 


2 


2/im" 3 


Critical concentration, c cr 


1 


1 /im -3 


Effective diffusion constant, D = Ma 


0.007 


0.7/im 2 /s 


Rotational viscosity, T 


1 


1 kPa • s 


Friction coefficient (if focal adhesion is present), 7 


- 0.04 


(0-0.04) nN-s//mi 4 


Activity, £ 


0-0.01 


(0 - 1) nN • /im 



